###############################################################
# Can Social Pressure Foster Responsiveness: Field Experiment
# by Bryant J. Moy
# Last Updated: 5/28/20 
###############################################################

rm(list=ls()) # Clear local environment
library(tidyverse);library(sandwich);library(AER);library("survminer");
library(knitr);library(broom);library(pwr);library(aod);library(dotwhisker)

##################################################################
# This public data file does not include the control variables, ##
# because I made a commitment to the IRB to keep the public     ##
# officials in the experiment anonymous. Consequently, this R    ##
# file does not recreate the tables in the appendices.          ##
# Results from the paper can be replicated using lines 20 - 56. ##
# The manipulation check starts at line 350                     ##
##################################################################

load("Moy_JEPS.RData")

####################################################
# Table 2. Response Rates by Treatment Condition
####################################################
data %>%
  group_by(AssignedT) %>%
  summarise_at(vars(responsive),mean) 
# Control .535 
# Treat_duty .550
# Treat_Peer .467

data %>%
  group_by(AssignedT)%>% 
  filter(responsive==1)%>%
  summarize(n = n())
# Treat_Duty   258
# Control      251
# Treat_Peer   220

data %>%
  group_by(AssignedT)%>% 
  summarize(n = n())
# Treat_Duty   469
# Control      469
# Treat_Peer   471


##########
# Table 3
##########
# ITT
ITT<-lm(responsive~Peer+Duty,data=data)
summary(ITT,vcovHC(ITT),correlation = F)

# IV
cace <- ivreg(responsive~Peer_Open+Duty_Open,~Peer+Duty, data=data) # z is instruments
coeftest(cace, vcovHC(cace),correlation = F)


################################################
# Appendix
# See note starting on line 10.
# Jump to line 148 for the power analysis
# Jump to line 350 for the manipulation check 
###############################################

############
# Table A1
#############
vars <- c("Population", "City_Age", "City_Female","City_Black","City_White","City_Renter")

v_n<- lapply(data[vars], NROW)
v_mean <- lapply(data[vars], mean)
v_mean <- round(as.numeric(v_mean), 1)
v_sd <- lapply(data[vars], sd)
v_sd <- round(as.numeric(v_sd), 1)
v_min <- lapply(data[vars], min)
v_min <- round(as.numeric(v_min), 1)
v_max <- lapply(data[vars], max)
v_max <- round(as.numeric(v_max), 1)

v_tab <- cbind(N=v_n,mean=v_mean, St_Dev=v_sd, min=v_min, max=v_max)

kable(v_tab)
  

##########################
# Table A2 and Table A3
##########################

## Table A3
mylogit_Duty <- glm(Duty~Population + City_Age+City_Female+ City_Black+City_White+ City_Renter, data = data, family = "binomial")
mylogit_Peer <- glm(Peer~Population + City_Age+City_Female+ City_Black+City_White+City_Renter, data = data, family = "binomial")
mylogit_Control <- glm(Control~ Population + City_Age+City_Female+ City_Black+City_White+City_Renter, data = data, family = "binomial")
stargazer:: stargazer(mylogit_Duty, mylogit_Peer,mylogit_Control, type= "text")

## Table A2

#Overal effect of observables on the assignment are not statistically significant. P(> X2) are not significant
wald.test(b = coef(mylogit_Duty), Sigma = vcov(mylogit_Duty), Terms = 2:7) # .42
pchisq(6, df=6, lower.tail = F) #0.4231901
wald.test(b = coef(mylogit_Peer), Sigma = vcov(mylogit_Peer), Terms = 2:7) # .58
pchisq(4.8, df=6, lower.tail = F) #0.5697087
wald.test(b = coef(mylogit_Control), Sigma = vcov(mylogit_Control), Terms = 2:7) #.86
pchisq(2.6, df=6, lower.tail = F) #0.8571125

Table_A2<- tribble(
  ~treatment, ~prob, ~n,
  "Duty", 0.423, 469 ,
  "Peer", 0.570, 471,
  "Control", 0.8571125, 469
)

Table_A2


###########
# Table B1
###########
data %>%
  filter(OPEN==1) %>%
  group_by(AssignedT)%>% 
  summarize(n = n())
# Treat_Duty   376
# Control      364
# Treat_Peer   367

data %>%
  group_by(AssignedT)%>% 
  summarize(n = n())
# Treat_Duty   469
# Control      469
# Treat_Peer   471

data %>%
  group_by(AssignedT)%>% 
  summarise_at(vars(OPEN),mean)
# Treat_Duty   77.6
# Control      80.2
# Treat_Peer   77.9

############
# Table B2
############
conditioned<-lm(responsive~Peer+Duty+OPEN,data=data)
summary(conditioned,vcovHC(conditioned),correlation = F)


###################
# Power Analysis
###################

C_1<-pwr.anova.test(k = 3, n = NULL, f = .1, sig.level = 0.05, power = .8)
C_1
plot(C_1)


C_2<-pwr.anova.test(k = 3, n = 467, f = NULL, sig.level = 0.05, power = .8)
C_2
plot(C_2)



#################################
# For Duration1: initial
# Risk ends at day 68
# after examining table(data$duration1), the last day is 68 day.
# I code all NA (which stand for non-response or not opened email) to 70 (68+2)
####################################
table(data$duration1)

# Duration
qplot(data$duration1,geom="bar", main = "Days to Initial Response",#fill=I("grey"),
      ylab = "Count of Mayoral Responses", xlab = "Days")+ xlim(c(0,70))+ylim(c(0,140))+theme(plot.title = element_text(hjust = 0.5))
# Note: Not shown are the 681 coded as non-response after 68 days

# Analysis of Duration1: Initial Response
fit<-survfit(Surv(duration1, responsive) ~ AssignedT, data = data)
res.cox <- coxph(Surv(duration1, responsive) ~ AssignedT, data = data)
res.cox
# 17% low (1-.83) 
# Relative to baseline, the hazard of response is 17% lower in the peer condition rather than control
stargazer::stargazer(res.cox,type="text")

ggsurvplot(fit, conf.int = F,  title="Survival Curve for Initial Responses",
           risk.table = TRUE, xlim= c(0,72),
           palette= "grey", legend.labs=c("Control", "Duty", "Peer"),
           ggtheme = theme_minimal()+theme(plot.title = element_text(hjust = 0.5)), data = data)


################################
# Heterogenous Effects
################################

# Variables are collected from different places like the Better Government Association and the National Freedom of Information Coalition. 
# As such, there is wide missingness in covariates outside of city demographics listed in Table A1. The following is for exploration purposes, 
# not causal inference.

# Peer
Peer_strongT<-tidy(lm(responsive~log(Population)+
             City_Age+City_Female+City_Black+
             City_Renter+Peer*strong_time1+
            +cp_fog+cp_partisan+
              cp_termlength+cp_termlimits,
           data=data[which(data$Peer==1 | data$Control==1),]))[15,]

Peer_strongM<-tidy(lm(responsive~log(Population)+
             City_Age+City_Female+City_Black+
             City_Renter+strong_time1+
             +Peer*cp_fog+cp_partisan+
             cp_termlength+cp_termlimits,
           data=data[which(data$Peer==1 | data$Control==1),]))[15,]

Peer_Partisan<-tidy(lm(responsive~log(Population)+
             City_Age+City_Female+City_Black+
             City_Renter+strong_time1+
             +cp_fog+Peer*cp_partisan+
             cp_termlength+cp_termlimits,
           data=data[which(data$Peer==1 | data$Control==1),]))[15,]

Peer_termlimit<-tidy(lm(responsive~log(Population)+
             City_Age+City_Female+City_Black+
             City_Renter+strong_time1+
             +cp_fog+cp_partisan+
             cp_termlength+Peer*cp_termlimits,
           data=data[which(data$Peer==1 | data$Control==1),]))[15,]

Peer_termlength<-tidy(lm(responsive~log(Population)+
             City_Age+City_Female+City_Black+
             City_Renter+strong_time1+
             +cp_fog+cp_partisan+
             Peer*cp_termlength+cp_termlimits,
           data=data[which(data$Peer==1 | data$Control==1),]))[15:17,]


Peer_Party<-tidy(lm(responsive~log(Population)+
             City_Age+City_Female+City_Black+
             City_Renter+strong_time1+
             cp_fog+cp_partisan+Peer*Party_Democrat+Peer*Party_Republican+
             cp_termlength+cp_termlimits,
           data=data[which(data$Peer==1 | data$Control==1),]))[17:18,]

Peer_Het<-rbind(Peer_strongT,Peer_strongM,Peer_Partisan,Peer_termlimit,Peer_termlength,Peer_Party)

Peer_Het$term

# rename the terms
Peer_Het$term<-c("Strong Time", "Strong Mayor", "Partisan", "Term Limit","Term Length: 2 yrs","Term Length: 3 yrs", 
                          "Term Length: 4 yrs","Democrat",
                          "Republican")

dwplot(Peer_Het, 
        vline = geom_vline(xintercept = 0, colour = "grey60", linetype = 2)) + # plot line at zero _behind_ coefs
   theme_grey() + xlab("Treatment Interaction Effects") + ylab("") +
   #ggtitle("Interaction of Peer Effects on Responsiveness") +
   theme(plot.title = element_text(face="bold"), legend.position = "none")



# Duty
Duty_strongT<-tidy(lm(responsive~log(Population)+
                        City_Age+City_Female+City_Black+
                        City_Renter+Duty*strong_time1+
                        +cp_fog+cp_partisan+
                        cp_termlength+cp_termlimits,
                      data=data[which(data$Duty==1 | data$Control==1),]))[15,]

Duty_strongM<-tidy(lm(responsive~log(Population)+
                        City_Age+City_Female+City_Black+
                        City_Renter+strong_time1+
                        +Duty*cp_fog+cp_partisan+
                        cp_termlength+cp_termlimits,
                      data=data[which(data$Duty==1 | data$Control==1),]))[15,]


Duty_Partisan<-tidy(lm(responsive~log(Population)+
                         City_Age+City_Female+City_Black+
                         City_Renter+strong_time1+
                         +cp_fog+Duty*cp_partisan+
                         cp_termlength+cp_termlimits,
                       data=data[which(data$Duty==1 | data$Control==1),]))[15,]

Duty_termlimit<-tidy(lm(responsive~log(Population)+
                          City_Age+City_Female+City_Black+
                          City_Renter+strong_time1+
                          +cp_fog+cp_partisan+
                          cp_termlength+Duty*cp_termlimits,
                        data=data[which(data$Duty==1 | data$Control==1),]))[15,]


Duty_termlength<-tidy(lm(responsive~log(Population)+
                           City_Age+City_Female+City_Black+
                           City_Renter+strong_time1+
                           +cp_fog+cp_partisan+
                           Duty*cp_termlength+cp_termlimits,
                         data=data[which(data$Duty==1 | data$Control==1),]))[15:17,]


Duty_Party<-tidy(lm(responsive~log(Population)+
                      City_Age+City_Female+City_Black+
                      City_Renter+strong_time1+
                      cp_fog+cp_partisan+Duty*Party_Democrat+Duty*Party_Republican+
                      cp_termlength+cp_termlimits,
                    data=data[which(data$Duty==1 | data$Control==1),]))[17:18,]


Duty_Het<-rbind(Duty_strongT,Duty_strongM,Duty_Partisan,Duty_termlimit,Duty_termlength,Duty_Party)

Duty_Het$term

Duty_Het$term<-c("Strong Time", "Strong Mayor", "Partisan", "Term Limit","Term Length: 2 yrs","Term Length: 3 yrs", 
           "Term Length: 4 yrs","Democrat",
           "Republican")

dwplot(Duty_Het, 
       vline = geom_vline(xintercept = 0, colour = "grey60", linetype = 2)) + # plot line at zero _behind_ coefs
  theme_grey() + xlab("Interaction of Duty on Responsiveness") + ylab("") +
  #ggtitle("Interaction of Social Accountability on Responsiveness") +
  theme(plot.title = element_text(face="bold"), legend.position = "none")



#########################
# P-Value Corrections
#########################
### P-Value Adjust
pvals_ITT<-c(summary(ITT,vcovHC(ITT), correlation = F)$coefficients[2:3,4])
pvals_CACE<-c(summary(cace,vcovHC(cace), correlation = F)$coefficients[2:3,4])
BONF_ITT = p.adjust(pvals_ITT, "bonferroni")
BH_ITT = p.adjust(pvals_ITT, "BH")
ITT_res = cbind(pvals_ITT, Bonferroni=BONF_ITT,BH=BH_ITT)
ITT_res
pvals_CACE<-c(summary(cace)$coefficients[2:3,4])
BONF_CACE = p.adjust(pvals_CACE, "bonferroni")
BH_CACE = p.adjust(pvals_CACE, "BH")
CACE_res = cbind(pvals_CACE, Bonferroni=BONF_CACE,BH=BH_CACE)
CACE_res
# In Table
tibble(Model = c("ITT","ITT","CACE","CACE"), Treatment= c("Duty","Peer","Duty","Peer"), 
       "Raw P-Value" = c(ITT_res[2,1], ITT_res[1,1],CACE_res[2,1], CACE_res[1,1]), 
       Bonferroni = c(ITT_res[2,2],ITT_res[1,2],CACE_res[2,2],CACE_res[1,2]), 
       BH =c(ITT_res[2,3],ITT_res[1,3],CACE_res[2,3],CACE_res[1,3]))







################################
# Manipulation Check: Survey
################################

library(likert)
require(grid)
require(lattice)
require(latticeExtra)

MC<- read_csv("Manipulation_Check_data.csv")
MC<-MC[-c(1,2),]


# Figure H1
MC$D1<-factor(MC$D1,
              levels = c("Not likely", "Somewhat likely", "Very likely"),
              ordered = T)
MC$P1<-factor(MC$P1,
              levels = c("Not likely", "Somewhat likely", "Very likely"),
              ordered = T)

h1<-data.frame(MC[,c(4,10)])
names(h1)<-c("Duty","Peer")
dimnames(h1)
results_h1<-likert(h1)
plot(results_h1,centered=F,include.center=T)



# Figure H2
MC$D3<-factor(MC$D3,
              levels = c("Strongly disagree", "Somewhat disagree", "Neither agree nor disagree","Somewhat agree","Strongly agree"),
              ordered = T)
MC$D4<-factor(MC$D4,
              levels = c("Strongly disagree", "Somewhat disagree", "Neither agree nor disagree","Somewhat agree","Strongly agree"),
              ordered = T)
MC$D5<-factor(MC$D5,
              levels = c("Strongly disagree", "Somewhat disagree", "Neither agree nor disagree","Somewhat agree","Strongly agree"),
              ordered = T)
h2<-data.frame(MC[,c("D3","D4","D5")])
names(h2)<-c("Non-Constituent: Less","Public's Belief: More","Overall Pressure")
results_h2<-likert(h2)
plot(results_h2,centered=T,include.center=T)


# Figure H3

MC$P3<-factor(MC$P3,
              levels = c("Strongly disagree", "Somewhat disagree", "Neither agree nor disagree","Somewhat agree","Strongly agree"),
              ordered = T)
MC$P4<-factor(MC$P4,
              levels = c("Strongly disagree", "Somewhat disagree", "Neither agree nor disagree","Somewhat agree","Strongly agree"),
              ordered = T)
MC$P5<-factor(MC$P5,
              levels = c("Strongly disagree", "Somewhat disagree", "Neither agree nor disagree","Somewhat agree","Strongly agree"),
              ordered = T)
MC$P6<-factor(MC$P6,
              levels = c("Strongly disagree", "Somewhat disagree", "Neither agree nor disagree","Somewhat agree","Strongly agree"),
              ordered = T)
MC$P7<-factor(MC$P7,
              levels = c("Strongly disagree", "Somewhat disagree", "Neither agree nor disagree","Somewhat agree","Strongly agree"),
              ordered = T)
MC$P8<-factor(MC$P8,
              levels = c("Strongly disagree", "Somewhat disagree", "Neither agree nor disagree","Somewhat agree","Strongly agree"),
              ordered = T)
h3<-data.frame(MC[,c("P3","P4","P5","P6","P7","P8")])
names(h3)<-c("Non-Constituent: Less","1400 Peers: Less","Report: More","Report to Peers: More","Email Access, Not Researcher","Overall Pressure")

results_h3<-likert(h3)
plot(results_h3, centered=T,include.center=T) 



